
#ifndef SERIESEPSILON_H
#define SERIESEPSILON_H

#include "SeriesBase.h"

// A series extrapolator using Wynn's Epsilon Algorithm

template <class T> 
class SeriesEpsilon : public virtual SeriesBase<T>
{
#define INITVECS(n) \
    do { \
        s.reserve(n); \
        t.reserve(n); \
        if (doerrcalc) { \
            ds.reserve(n); \
            dt.reserve(n);  \
        } \
    } while (false)

public:
    SeriesEpsilon(uint n=20, bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<T>(n,_doerrcalc,_saveterms), i0(0), startsize(n) 
    { INITVECS(n); }

    SeriesEpsilon(const vector<T>& v, bool _doerrcalc=true,
                  bool _saveterms=true) : 
        SeriesBase<T>(v.size(),_doerrcalc,_saveterms), i0(0),
        startsize(v.size()) 
    { INITVECS(v.size()); AddTerms(v); }

    SeriesEpsilon(const vector<T>& v, const vector<T>& verr, 
                  bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<T>(v.size(),_doerrcalc,_saveterms), i0(0),
        startsize(v.size()) 
    { INITVECS(v.size()); AddTerms(v,verr); }

    ~SeriesEpsilon() {}
#undef INITVECS

protected:
    void CalcResult(const T& term);
    const T& GetResult() const {return s[i0];}
    const vector<T>& GetDSum() const {return ds[i0];}

private:
    vector<T> s;  // A series of approximations to the sum
    // The latest best guess is s[i0]
    vector<T> t;  // A list of helper quantities
    vector<vector<T> > ds; // ds[i][j] = d(s[i])/d(term_j)
    vector<vector<T> > dt; // dt[i][j] = d(t[i])/d(term_j)
    uint i0;
    uint startsize;
};

template <class T>
inline void SeriesEpsilon<T>::CalcResult(const T& term)
{
    uint n=NTerms();
    assert(s.size()+1 == n);
    assert(n==1 || t.size()+2 == n);
    if (doerrcalc) {
        assert(ds.size() == s.size());
        assert(dt.size() == t.size());
    }

    // Update the tables for the new term
    if (n == 1) {
        s.push_back(term);
        if (doerrcalc) {
            ds.push_back(vector<T>(1,T(1.0)));
            ds.back().reserve(startsize);
        }
        i0 = 0;
    } else {
        s.push_back(PlainSum());
        T invterm = 1.0/term;
        t.push_back(invterm);
        if (doerrcalc) {
            for (uint i = 0; i < n-2; i++) {
                ds[i].push_back(T(0.0));
                dt[i].push_back(T(0.0));
            }
            ds[n-2].push_back(T(0.0));
            ds.push_back(vector<T>(n,1.));
            ds.back().reserve(startsize);
            dt.push_back(vector<T>(n,0.));
            dt.back().reserve(startsize);
            dt.back()[n-1] = -invterm*invterm;
        }

        int i;
        for (i = n-3; i >= 0; i-=2) {
            T invtdiff = 1.0/(t[i+1]-t[i]);
            s[i] = s[i+1] + invtdiff;
            const T invtdiffsq = invtdiff*invtdiff;
            if (doerrcalc) {
                for (int j = 0; j < n; j++) {
                    ds[i][j] = ds[i+1][j] - invtdiffsq*(dt[i+1][j]-dt[i][j]);
                }
            }
            if (i>0) {
                T invsdiff = 1.0/(s[i]-s[i-1]);
                t[i-1] = t[i] + invsdiff;
                if (doerrcalc) {
                    for (uint j=0; j<n;j++) {
                        dt[i-1][j] = dt[i][j] - invtdiffsq*(ds[i][j]-ds[i-1][j]);
                    }
                }
            }
        }
        i+=2;
        assert(i==0 || i==1);
        i0 = i;
    }
}

#endif
